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Recently B. Derrida [Phys. Rev. E 55, 3705 (1997)] introduced a numerical technique that 
allows one to measure the fraction of persistent spins in a coarsening nonequilibrium system at 
finite temperature. In the present work we extend this method in a way that domain walls can 
be clearly identified. To this end we consider three replicas instead of two. As an application 
we measure the surface area of coarsening domains in the two-dimensional Ising model at finite 
temperature. We also discuss the question of to what extent the results depend on the algorithmic 
implementation. 

I. INTRODUCTION 

Dynamical systems quenched from a disordered into an ordered phase may display interesting coarsening phenom- 
ena A simple example is the Ising model evolving by heat bath (HBD) or Glauber dynamics (GD). In the ordering 
phase of this model, starting with random initial conditions, patterns of ferromagnetic domains are formed whose 
typical size grows with time as t 1 / 2 . For zero temperature, these domains are fully ordered and the domain walls 
evolving in time can be identified as bonds between oppositely oriented spins ||. For nonzero temperature, however, 
it is hard to define domains and domain walls because it is difficult to distinguish between 'true' domains and minority 
islands generated by thermal fluctuations. This situation emerges, for example, in the two-dimensional (2D) Ising 
model at finite temperature below T c (see Fig. |l|a). 

Recently B. Derrida || proposed a method that allows one to measure properties related to coarsening in presence 
of thermal fluctuations. The main idea of this method lies in the comparision of two identical copies (replicas) A 
and B of the same system. Both replicas are submitted to the same thermal noise, i.e., their numerical updates 
are determined by the same sequence of random numbers. Copy A starts with random initial conditions and begins 
to coarsen whereas copy B starts from a fully magnetized state and therefore remains ordered as time evolves. The 
assumption is that all spin flips occurring in replica B can be regarded as thermal fluctuations. Therefore, when a spin 
flip occurs simultaneously in both replicas, it can be considered as a thermal fluctuation, otherwise as a fluctuation 
due to the coarsening process. 

In Refs. HQ this method was used to determine the fraction of persistent spins Q at nonzero temperature as a 
function of time. At zero temperature a spin is said to be persistent up to time t if it never flipped before. In the 
Ising model the fraction of persistent spins r(t) decays according to a power law r(t) ~ t~ e where 9 is an independent 
exponent. For the Id Glauber model it was proved that 9 = 3/8 || whereas in higher dimensions 9 could only be 
determined by numerical simulations and approximation methods H. For T > 0, however, the fraction of spins 
that never flipped decays exponentially since thermal fluctuations occur everywhere at some finite rate. To overcome 
this difficulty, Derrida proposed to consider a spin as "persistent" if its temporal evolution in copies A and B is fully 
synchronized up to time t. Using this definition of persistence he analyzed the 2D Ising model and observed that r(t) 
decays algebraically for < T < T c and saturates at some finite value for T > T c . Below the critical temperature the 
exponent 9 seems to be the same as for T = while at T = T c a different exponent is observed. 

An imperfection of the method developed by Derrida is that only one type of domain can be identified, namely, 
those that are magnetized in the same way as system B. Therefore different spin flips in A and B indicate the presence 
of oppositely magnetized domains rather than the presence of a domain wall. This means that persistent spins can be 
identified only in those domains which have the same orientation as copy A. For the same reason the method cannot 
be used to analyze other properties such as, for example, the dynamics of domain walls. 

In the present work we extend Derrida's method in a way that domain walls can be identified. For this purpose 
we consider three replicas A, B, C instead of two. As before, all replicas are submitted to the same realization of 
noise. Replica A starts with random initial conditions and serves as the master copy in which the coarsening process 
takes place. The temporal evolution of replica A is compared with that of replicas B and C, which start from fully 
ordered initial conditions with positive and negative magnetization, respectively. As in the original setup, domains 
with positive magnetization in copy A exhibit the same thermal fluctuations as copy B. Likewise, thermal fluctuations 
in domains with negative magnetization in copy A are synchronized with those in copy C . Along the domain walls, 
however, fluctuations in replica A may occur that are different from those in B as well as in C. Detecting such 
fluctuations by an appropriate observable (to be defined below) we are able to identify domain walls in a coarsening 
process at nonzero temperature. The remarkable efficiency of this method is illustrated in Fig. [lib. In addition, our 
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FIG. 1. The 2D Ising model evolving by heat bath dynamics, (a): of a system with 250 x 250 sites and periodic boundary 
conditions at temperature T = 0.9 T c after 500 updates, (b): Corresponding domain walls detected by the observable A; 
defined in Eq. ([j]). The arrows point at missing links of the contours. 



technique allows one to measure interesting physical quantities such as, for example, the surface area of domains as a 
function of time. In what follows we restrict ourselves to the 2D Ising model evolving by HBD and GD. However, the 
technique can easily be generalized and may be applied to many other stochastic coarsening processes. For example, 
applying the method to the Potts model with q > 2 states per site requires introducing q + 1 different replicas. 

A fundamental problem of numerical methods based on several replicas evolving under the same noise is that the 
results may depend on the algorithmic implementation. This was first observed in so-called damage spreading (DS) 
problems. In DS simulations J§| two replicas of a nonequilibrium system, submitted to the same thermal noise, are 
started from slightly different initial conditions. If the difference between the two copies (the damage) stays finite 
or even diverges, the system is said to exhibit damage spreading. Otherwise, if the two replicas merge into a fully 
synchronized evolution, damage is said to heal. Initially DS fascinated researchers, since it would have indicated the 
existence of different dynamical phases in stochastic models analogous to chaotic and regular phases in deterministic 
systems. However, later it was realized Q that such DS phases are ambiguous since the usage of different but 
equivalent algorithms for the same dynamical system can lead to different DS phase structures [jlOfl - For example, 
in the Ising model with HBD damage always heals while in the case of GD damage may spread |n[. The reason is 
that GD and HBD, although indistinguishable on a single replica, are characterized by different correlations when 
two or more replicas are simulated using the same random numbers [Q. As we are going to demonstrate, a similar 
algorithmic dependence appears in the present numerical technique where several replicas submitted to the same noise 
are used to analyze coarsening processes. Thus we have to verify to what extent the results obtained by Derrida || 
in the case of HBD [j] are physically relevant or rather artifacts of different algorithmic schemes. 

The article is organized in the following way. In Sec. H we define HBD and GD as well as an observable by which 
domain walls can be detected. Our numerical results for HBD are presented in Sec. III. By comparing results for 
HBD and GD, we address the problem of algorithmic independence in Sec. IV. Finally our results are summarized 
and discussed in Sec. |v|. 



II. DETECTION OF DOMAIN WALLS IN THE ISING MODEL 

The Ising model evolving by HBD or GD is defined as follows. Consider a (/-dimensional square lattice with spins 
Oi(t) = ±1. The energy at time t is given by 



1 According to the usual terminology the dynamical rule used by Derrida in Ref. 
dynamics rather than Glauber (spin flip) dynamics. 



is denoted as heat bath (spin orienting) 
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where j runs over the nearest neighbors of site i. The local field hi(t) determines the transition probability Pi(t) for 
the spin a% at time t: 

e hi{t)/k B T 

= e hi{t)/k B T + e -hi(t)/k B T- ( 2 ) 

HBD and GD differ in their update rules: In HBD the spin <7j is oriented according to the local field hiit) by 

a i {t + l)=sgn\p i {t)-z i {t)], (3) 

where Zi(t) are independent random numbers drawn from a uniform distribution between and 1. On the other hand, 
in GD the spin <7j is flipped depending on its previous orientation: 

, . ^ _ / o-i(t)sgn\pi(t) - Zi(t)] if Ui{t) = +1 , , 

^ + iJ ~ 1 ^sgnri-^Ci) - *(t)] if ^(t) = -1 ' ^ 

One can easily verify that in both dynamics the probability to get <7j(f + 1) = +1 is the same, as expected from the 
equivalence between HBD and GD. 

Let us now consider three replicas A, B, and C and denote their spins by crf(t), crf(t), and erf (t). As stated before, 
the initial conditions are 

^(0)=sgn[i-4 0) ], of(0)=+l, <7p(0) = -l, (5) 

where z^ are random numbers between and f . The three replicas evolve under the same realization of noise, i.e., 
the same random numbers Zi(t) are used for the updates of crf(t), crf(t), and <jf(t). 

We now turn to the observables we want to analyze. Derrida's definition for the fraction of persistent spins r(t) 
can be generalized easily: A spin af(t) is said to be "persistent" up to time t if it experienced exclusively thermal 
fluctuations, which means that it was synchronized either with o~f(t) or with o~f(t) for the whole time, i.e., 

l + af{t)af{t) , ^ l + a?{t)af{t) 




n * o - • w 



0<t'<t 0<t'<t 



where N is the total number of sites. We will analyze this quantity numerically in Sees, pi -IV. 

In order to identify domain walls, we define an observable Aj(i) which compares replicas A, B, and C at site i and 
its nearest neighbors. We consider site i as belonging to a domain wall (i.e., Aj(i) = 1) if (a) site i or at least one of 
its nearest neighbors in copies A and B are in different states, and (b) if site i or at least one of its nearest neighbors 
in copies A and C are in different states. Since HBD and GD evolve independently on two (even and odd) sublattices, 
we assume that these nearest neighbors belong to the same sublattice (for example, if i — (x, y), the nearest neighbors 
on the same sublattice are (x ± 2, y) and (x, y ± 2)). Formally the observable Ai(t) is defined by 

where j runs over site i and its nearest neighbors on the same sublattice. It turns out that this observable allows one 
to identify domain walls, as illustrated for HBD in Fig. |l|b. 

The definition (|) appears to be quite complicated since it involves the nearest neighbors of site i. It would have 
been more natural to define a local observable 5i(t), which is 1 if o~f is different from of and o~f (indicating a 
fluctuation generated by the coarsening process) and otherwise: 

m= (i^imM.)(i^iMB). (8) 



But, using the initial conditions specified in Eq. (||), it would turn out that Si(t) = for all t and for both HBD and 
GD. This is due to an overlap of the regions in A which are synchronized with either B or C. For HBD this can be 
proven as follows. Assume that the three replicas at time t are in a state where the inequality 
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FIG. 2. Heat bath dynamics. Measurement of the fraction of persistent spins r(t) at various temperatures as a function 
of time (c.f. Ref. [3]). The bold line indicates the slope —0.22. The right hand graph shows a data collapse of the rescaled 
quantity r(t)/e 2 in the supercritical regime T > T c . 



P?(t)<P?(t)< P ?(t) (9) 

holds for all i. Since the number of positive spins generated by the update rule (0) for a given random number Zi(t) is 
monotonically increasing with Pi(t), one can show that hf (t + 1) < hf(t + 1) < hf (t + 1). Therefore the inequality (§) 
is also satisfied at the next time step t + 1 . Since this inequality is satisfied by the initial conditions (|5|) , it holds by 
induction at any time. This implies that events with of ^ of and of ^ of do not occur, hence Si(t) = for HBD. 
In the case of GD the proof is trivial: Since replicas B and C evolve precisely in opposite states (of(t) = — erf (i)), 
the observable Si(t) vanishes automatically. Thus the local observable defined in Eq. (0) cannot be used in order to 
detect domain walls. This is the reason we use the more complicated definition of Eq. M). 



III. HEAT BATH DYNAMICS: PERSISTENT SPINS AND THE FRACTAL DIMENSION OF DOMAIN 

WALLS 

In this section we present numerical results for the 2D Ising model with HBD. We simulate three replicas of a 
system of size 1000 x 1000 with periodic boundary conditions. Starting with the initial conditions ((HJ) we measure the 
fraction of persistent spins r(t) as defined in Eq. (||) and the average Peierls length (circumference) of the domains 

a(*) = ^ ( 10 ) 

i 

The quantities r(t) and A(t) are measured up to 5000 time steps and averaged over 10 independent runs. Our 
simulation data are shown in Figs. ^M. 

The results for the fraction of persistent spins r(t) (see left hand graph in Fig. ||) are in full agreement with Ref. [||. 
For T = wc observe an algebraic decay r(t) ~ t~ e with 9 = 0.23 ± 0.03. For < T < T c r(t) decreases rapidly on 
short time ranges while for longer times it crosses over to a T-independent power law decay. Precisely at the critical 
temperature, however, r(t) seems to vanish like r(t) ~ t~ e " with a different exponent 9 C ~ 1 (the physical relevance 
of this exponent will be discussed in the Sec. ^vj ). Finally, for T > T c , r(t) saturates at some finite value. This can 
be explained as follows. For T > T c the total magnetization in copies B and C decays exponentially. It has been 
shown fi"l|| that under these conditions any difference between two replicas evolving by HBD vanishes exponentially, 
i.e., damage heals spontaneously. This means that the replicas A,B,C converge and eventually merge into a fully 
synchronized evolution within finite time and consequently a finite fraction of persistent spins survives. Notice that 
in the limit T — > oo all replicas are already synchronized after a single time step. For finite temperatures we observe 
a scaling behavior (see right hand graph in Fig. ||) 

r(t) ~ (T-T c ) z f(t(T-T c ) z ), (T>T C ) (11) 

where z ~ 2.125 is the dynamical critical exponent of HB dynamics and fix) is a scaling function which behaves as 
f(x) ~ 1/x for x — > and saturates for x — > oo. 
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The results for the total Peierls length of the domains A(t) illustrated in Fig. [3| indicate a power law behavior 
A(t) = with k ~ 0.5 in the regime < T < T c and an exponential decay in the disordered phase T > T c (to 
determine the behavior at T = T c more numerical effort would be needed, c.f. Ref. that can be explained by the 
synchronization of the copies within finite time. The corresponding scaling behavior reads 

A(t) ~ (T - T C ) KZ g(t(T - T c ) z ) , (T > T c ) (12) 

where g(x) is a scaling function as shown in the right hand graph in Fig. ^. 

The exponent k ~ 0.5 in the subcritical regime may be interpreted as follows. The total number of domains n(t) 
in a large but finite sized system decreases as n(t) ~ t~ d l 2 . Moreover, the average size of the domains £ grows with 
time as £ ~ t 1 ^ 2 . Hence in the 2D Ising model the Peierls length behaves as 

n(t)£ ~ A(t) ~ r 1 ' 2 . (13) 

This result suggests that surfaces of domains in a coarsening Ising model are regular, i.e., they do not have a fractal 
structure at finite temperatures T < T c . This result can be explained as follows. The coarsening process is driven by 
the tendency of the systen to minimize its energy H, i.e., to minimize the surface area (the Peierls length A(£)) of 
the domains. This makes it highly unlikely for the domain walls to form fractal structures. This mechanism works 
not only at zero temperature but prevails in the entire subcritical regime. To support this argument we plotted the 
energy H against the Peierls length in Fig. |^. For T <T C the curves are monotonically decreasing with time and seem 
to have a well defined minimum at A = 0. This suggests that the mechanism for coarsening, apart from different time 
scales, is the same in the whole subcritical regime. In the disordered phase, however, the curves have a flat shape 
close to A = and therefore the dynamics of A(t) is no longer driven by the minimization of energy. 



IV. THE PROBLEM OF ALGORITHMIC DEPENDENCE 



As outlined in the Introduction, any numerical technique based on several replicas of a nonequilibrium system 
evolving under the same realization of noise may depend in a crucial way on the algorithmic implementation of the 
dynamics. We now discuss this dependence in the present problem by comparing HBD and GD. 

HBD and GD are two different but equally legitimate algorithmic implementations of the same nonequilibrium 
process that mimics the evolution of an Ising model in contact with a thermal reservoir. To understand this, it may 
be helpful to rewrite the update rule for GD (Q) as 

n n , n _ / sgn[pi(f) - Zi(t)] ifoi(t) = +l , , 

H + j I sgn[ Pl (t) - {1 - Zi (t)}} if <Ti(t) = -1 ' li4j 

This rule differs from HBD only inasmuch - depending on ai(t) - the random number Zi(t) or 1 — Zi(t) is used. Since 
in a simulation of a single replica each random number is used only once, it makes no difference whether Zi(t) or 
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FIG. 4. Total energy H(t) of a coarsening Ising model as a function of the Peierls length A(t) for various temperatures. 



1 — Zi(t) enters the update rule and therefore both dynamical procedures are fully equivalent. In other words, looking 
at the space-time trajectory of a single Ising model, one cannot distinguish whether it was generated by HBD or 
GD. However, if we consider more than one replica evolving under the same noise, each random number Zi(t) is used 
several times leading to different correlations between the replicas depending on whether HBD or GD are used. 

HBD and GD are two members of an infinite family of equivalent dynamical rules jl3|,[l2| . All these rules are equally 
legitimate and there is no reason to prefer particular rules such as HBD or GD. Physical properties, however, should 
not depend on the choice of the dynamical procedure. Thus, in order to prove the algorithmic independence of a 
specific result, one would have to verify all these rules separately. Since this is practically impossible, we will restrict 
ourselves to the examples of HBD and GD. We will show that some of the previous results are affected by a change 
of the procedure and hence are physically irrelevant whereas others are not. 

We now repeat the numerical simulations described in the previous section using GD instead of HBD. A snapshot 
of the simulation is shown in Fig. ^[ Comparing Figs, ^b andj^b we notice that the HBD algorithm yields very clean 
shapes for the Peierls contours while GD produces a lot of additional fluctuations. On the other hand the contours 
in the Glauber case are strictly closed because of symmetry reasons whereas for HBD there are missing links (two of 
them are marked by arrows in Fig. |l|b). 

We would like to note that any dynamical rule used in the present problem should synchronize the evolution of 
replicas A, B or A, C in the interior of large domains. In the terminology of DS this means that damage has to 
heal. As HBD is the most correlated algorithm, damage heals very rapidly and thus yields high resolution in the 
determination of the walls. On the other hand GD in 2D is known to exhibit a DS transition at T = T s ~ 0.95 T c p4[ . 
The simulations in Fig. at T = 0.9 T c take place very close to this transition which explains why the output is 
rather noisy. 

The numerical results for GD and HBD arc significantly different (see Fig. ^J). Only at zero temperature where 
GD and HBD coincide does one obtain identical results. For < T < T c the fraction of persistent spins r(t) exhibits 
similar behavior in both cases: It first decays and then crosses over to an algebraic decay r(t) ~ t~ e with the same 
~ 0.2 as in the T = case. This suggests that the crossover phenomenon is an intrinsic physical property of the 2D 
Ising model rather than an algorithmic artifact. At criticality, however, GD indicates an algebraic decay r(t) ~ i~ 2 . 
This is clearly different from the decay r(t) ~ t~ Y observed in the case of HBD. Therefore one may wonder about 
the physical meaning of r(t) at criticality. This problem may be circumvented by a systematic study of the crossover 
time as a function of temperature. 

The results for the total Peierls length of the domains are even more contradictory. While in the HBD case we 
found an algebraic decay A(t) ~ i -1 / 2 for < T < T e , we observe continuously varying exponents in the case of GD. 
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FIG. 5. The 2D Ising model evolving by Glauber dynamics, (a): Snapshot of a system with 250 x 250 sites and periodic 
boundary conditions at temperature T = 0.9 T c after 500 time steps, (b): Corresponding domain walls detected by the 
observable A;. Notice that all contours are strictly closed. 




If this was true, it would imply that the Peierls length per domain grows like A(t)/n(t) ~ £ K with some temperature- 
dependent exponent k > 1 indicating a "fractal" structure of the domain walls. Indeed, visual inspection of Fig. ^3 
makes it intuitively clear how such a "fractal" structure emerges. It is therefore tempting to discard the results for 
GD and to declare the smooth lines in Fig. [j]b as 'true' domain walls. However, as explained above, we have no 
justification for doing so! Therefore the method described in this paper cannot be used to determine the fractal 
dimension of domain walls on a safe ground. Nevertheless HBD gives us a lower bound for the exponent k. It is thus 
very likely, although not strictly proven, that domain walls in the coarsening Ising model at nonzero temperature do 
not have a fractal structure. 



V. CONCLUDING REMARKS 



By introducing three replicas of a 2D Ising model evolving under the same realization of noise we extend the 
numerical method proposed by B. Derrida || in a way that domain walls can be identified. Using HBD we measured 
the Peierls length of coarsening domains in the 2D Ising model. Our simulations confirm previous results and suggest 
that domain walls in the Ising model at T < T c are regular, i.e., they do not have fractal properties. 

A fundamental problem appearing here and related to the use of several replicas is the algorithmic dependence. As 
an example we compared HBD and GD. It turns out that the persistence exponent below the critical temperature is 
the same in both cases which suggests that this result is in fact related to a physical property rather than algorithmic 
artifacts. On the other hand, the Peierls length of the domain walls grows differently for HBD and GD, although it 
seems that HBD gives the correct result. 

In this context we would like to note that very recently a different method for the numerical estimation of the 
persistence exponent 6 has been proposed |l5| where the persistence probability is defined in terms of spin blocks. By 
analyzing the scaling behavior for different block sizes one can determine 9 even at finite temperature. Although this 
method cannot be used to identify domain walls, it is very interesting since it uses only a single replica wherefore the 
results do not depend on whether HBD or GD is used. In agreement with Ref. || and the present work the authors 
observe that the exponent 6 is the same in the entire subcritical regime. It would be interesting to use this method 
in order to determine 6 at criticality 
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